Supplmentary Material

Author

Co-author names

Published

January 15, 2024

Outstanding things and notes

1 Setting up

1.1 Loading packages

Code
pacman::p_load(tidyverse, # tidy family and related pacakges below
               kableExtra, # nice tables
               MuMIn,  # multi-model inference
               emmeans, # post-hoc tests
               gridExtra, # may not use this
               pander,   # nice tables
               metafor,  # package for meta-analysis
               ape,      # phylogenetic analysis
               MuMIn,  # multi-model inference
               patchwork,   # putting ggplots together - you need to install via devtool
               here,         # making reading files easy
               orchaRd # plotting orchard plots
               # more too add

)

eval(metafor:::.MuMIn)

1.2 Loading data: paper and tree data

Code
dat <- read.csv(here("data", "clean_data.csv"))

#dim(dat)
#head(dat)

tree <- read.tree(here("data", "species_tree.tre"))

1.3 Custom functions

Code
# function to calculate effect sizes
# Zr - correlation
# there is always n

effect_size <- function(m1, m2, sd1, sd2, n1, n2, n, # 12 arguments
                        est , se, p_val, direction, method){

  if(method == "mean_method"){
  
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s_pool <- sqrt( ((n1-1)*sd1^2 + (n2-1)*sd2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r    
    
  }else if(method == "count_method"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- sqrt(m1)
    s2 <- sqrt(m2)
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r     
    
  }else if(method == "percent_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "percent_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    sd1 <- sd1/100
    sd2 <- sd2/100
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1))
    m2 <- asin(sqrt(m2))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "t_method1"){
    

    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "t_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "F_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
  
  }else if(method == "F_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "correlation_method1"){
    
    r <- est
    
  }else if(method == "correlation_method2"){
    
    r <- 2*sin((pi/6)*est)
    
  }else if(method == "estimate_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  }else if(method == "estimate_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  } 
  
    if(r >= 1){
    # if over 1, we use 0.99 
    Zr <- atanh(0.99)
    }else if(r <= -1){
    Zr <- atanh(-0.99) # r = correlation
    } else {
    Zr <- atanh(r) # r = correlation
    } 
  
  VZr <- 1 /(n - 3)
  
  data.frame(ri = r, yi = Zr , vi = VZr)
  
}


##########
# functions for absolute values


# folded mean
folded_mu <-function(mean, variance){
  mu <- mean
  sigma <- sqrt(variance)
  fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma))
  fold_mu
} 

# folded variance
folded_v <-function(mean, variance){
  mu <- mean
  sigma <- sqrt(variance)
  fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma))
  fold_se <- sqrt(mu^2 + sigma^2 - fold_mu^2)
  # adding se to make bigger mean
  fold_v <-fold_se^2
  fold_v
} 



#' Title: Contrast name generator
#'
#' @param name: a vector of character strings
cont_gen <- function(name) {
  combination <- combn(name, 2)
  name_dat <- t(combination)
  names <- paste(name_dat[, 1], name_dat[, 2], sep = "-")
  return(names)
}

#' @title get_pred1: intercept-less model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred1 <- function (model, mod = " ") {
  name <- firstup(as.character(stringr::str_replace(row.names(model$beta), mod, "")))
  len <- length(name)
  
   if (len != 1) {
        newdata <- diag(len)
        pred <- metafor::predict.rma(model, 
                                     newmods = newdata,
                                     tau2.levels = 1:len)
    }
    else {
        pred <- metafor::predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title get_pred2: normal model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred2 <- function (model, mod = " ") {
  name <- as.factor(str_replace(row.names(model$beta), 
                                paste0("relevel", "\\(", mod,", ref = name","\\)"),""))
  len <- length(name)
  
  if(len != 1){
  newdata <- diag(len)
  pred <- predict.rma(model, intercept = FALSE, newmods = newdata[ ,-1])
  }
  else {
    pred <- predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title mr_results
#' @description Function to put results of meta-regression and its contrasts
#' @param res1: data frame 1
#' @param res1: data frame 2
mr_results <- function(res1, res2) {
  restuls <-tibble(
    `Fixed effect` = c(as.character(res1$name), cont_gen(res1$name)),
    Estimate = c(res1$estimate, res2$estimate),
    `Lower CI [0.025]` = c(res1$lowerCL, res2$lowerCL),
    `Upper CI  [0.975]` = c(res1$upperCL, res2$upperCL),
    `P value` = c(res1$pval, res2$pval),
    `Lower PI [0.025]` = c(res1$lowerPR, res2$lowerPR),
    `Upper PI  [0.975]` = c(res1$upperPR, res2$upperPR),
  )
}


#' @title all_models
#' @description Function to take all possible models and get their results
#' @param model: intercept-less model
#' @param mod: the name of a moderator 

all_models <- function(model, mod = " ", type = "normal") {
  
  # getting the level names out
  level_names <- levels(factor(model$data[[mod]]))
  dat2 <- model$data
  mod <- mod

  # creating a function to run models
  run_rma1 <- function(name) {
      
    # variance covarinace matrix for sampling errors
    VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, rho = 0.5)
    VCV1 <- nearPD(VCV1)$mat
    
      
    rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }

    run_rma2 <- function(name) {
    
            # variance covarinace matrix for sampling errors
            VCVa <- vcalc(vi = dat2$vi_abs, cluster = dat$shared_group, rho = 0.5)
            VCVa <- nearPD(VCVa)$mat
               
               rma.mv(abs_yi2, V = VCVa,
               mods = ~relevel(dat2[[mod]], ref = name),
                      random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
    }
    
    run_rma3 <- function(name) {
      
              # variance covarinace matrix for sampling errors
                VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, rho = 0.5)
                VCV1 <- nearPD(VCV1)$mat
               
                rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list( # hetero in relation to mod
                            formula(paste("~", mod, "| effectID")),
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     struct = "DIAG",
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }


# results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])
# this does not work for hetero model?
if (type == "normal"){

    model_all <- purrr::map(level_names[-length(level_names)], run_rma1)

  }else if(type == "absolute"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma2)
    
  }else if(type == "hetero"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma3)
    
  }
  
  # getting estimates from intercept-less models (means for all the groups)
  res1 <- get_pred1(model, mod = mod)
  
  # getting estiamtes from all contrast models
  res2_pre <- purrr::map(model_all, ~ get_pred2(.x, mod = mod))
  
  # a list of the numbers to take out unnecessary contrasts
  contra_list <- Map(seq, from=1, to=1:(length(level_names) - 1))
  res2 <- purrr::map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) 
  # creating a table
  res_tab <- mr_results(res1, res2) %>% 
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")
  
  # results
  res_tab

}

1.4 Preparing data: obtaining effect sizes etc.

Code
# calcuating effect sizes
effect2 <- pmap_dfr(list(dat$mean_group_1, dat$mean_group_2, 
                         dat$variance_group_1, dat$variance_group_2, 
                         dat$n_group_1, dat$n_group_2, dat$n, 
                         dat$effect_size_value, dat$effect_size_variance, 
                         dat$effect_size_p_value_numeric, 
                         dat$direction_change, dat$function_needed), 
                    effect_size)                    


# merging two data frames
dat <- cbind(dat, effect2)

# adding absolute value (effect sizes) too
dat <- dat %>% mutate(
  abs_yi = abs(yi), # we use this one (conservative)
  abs_yi2 = folded_mu(yi, vi), # alternative way
  abs_vi = folded_v(yi, vi)) %>% 
  mutate_if(is.character, as.factor) # changing all character to factor


# dat$Zr <- unlist(effect1[1,])
# dat$VZr <- unlist(effect1[2,])



# renaming X to effectID
colnames(dat)[colnames(dat) == "X"] <- "effectID"

# creating the phylogeny column

dat$phylogeny <-  gsub(" ", "_", dat$species_cleaned)

# checking species names match between tree and dataset

# match(unique(dat$phylogeny), tree$tip.label)
# match(tree$tip.label, unique(dat$phylogeny))
# 
# intersect(unique(dat$phylogeny), tree$tip.label)
# setdiff(unique(dat$phylogeny), tree$tip.label)
# 
# match(unique(dat$phylogeny), tree$tip.label)
# sum(is.na(match(unique(dat$phylogeny), tree$tip.label)))

# creating a correlation matrix from a tree


# creating a VCV for sampling error variances 
# we assure shared groups will have r = 0.5
VCV <- vcalc(vi = dat$vi, cluster = dat$shared_group, rho = 0.5)
VCV <- nearPD(VCV)$mat


# creating VCV for absolute values
VCV_abs <- vcalc(vi = dat$abs_vi, cluster = dat$shared_group, rho = 0.5)
VCV_abs <- nearPD(VCV_abs)$mat

2 Meta-analysis

3 Phylogenetic mulitlevle model

Code
# takes a while to run
mod <- rma.mv(yi = yi, 
              V = VCV,
              mod = ~ 1,
              data = dat,
              random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned,
                           ~ 1 | phylogeny),
              R= list(phylogeny = cor_tree),
              test = "t",
              sparse = TRUE)

# saving the runs
saveRDS(mod, here("Rdata", "mod.RDS"))
Code
# getting the save model
mod <- readRDS(here("Rdata", "mod.RDS"))

summary(mod)

Multivariate Meta-Analysis Model (k = 648; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-280.9830   561.9661   571.9661   594.3278   572.0597   

Variance Components:

            estim    sqrt  nlvls  fixed           factor    R 
sigma^2.1  0.1222  0.3495    648     no         effectID   no 
sigma^2.2  0.0028  0.0531    195     no          paperID   no 
sigma^2.3  0.0241  0.1553    139     no  species_cleaned   no 
sigma^2.4  0.0000  0.0000    139     no        phylogeny  yes 

Test for Heterogeneity:
Q(df = 647) = 387952763.7787, p-val < .0001

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0243  0.0292  -0.8306  647  0.4065  -0.0817  0.0331    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# hardly any phlogenetic effects so taking it out
round(i2_ml(mod),2)   
          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.64              80.00               1.85              15.80 
      I2_phylogeny 
              0.00 
Code
# phylogeney accounts nothing (0) so we can take it out
mod1 <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)

summary(mod1)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-263.5589   527.1177   535.1177   553.0195   535.1798   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1173  0.3425    650     no         effectID 
sigma^2.2  0.0051  0.0717    197     no          paperID 
sigma^2.3  0.0152  0.1232    140     no  species_cleaned 

Test for Heterogeneity:
Q(df = 649) = 366834657.1747, p-val < .0001

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0090  0.0273  -0.3305  649  0.7411  -0.0625  0.0445    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(i2_ml(mod1),2)   
          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.44              83.05               3.64              10.76 
Code
orchard_plot(mod1, xlab = "Effect Size: Zr", group = "paperID")

!!! - model comparison we should use ML - !!!

Code
# phylogeney accounts nothing (0) so we can take it out
mod1 <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)

summary(mod1)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-263.5589   527.1177   535.1177   553.0195   535.1798   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1173  0.3425    650     no         effectID 
sigma^2.2  0.0051  0.0717    197     no          paperID 
sigma^2.3  0.0152  0.1232    140     no  species_cleaned 

Test for Heterogeneity:
Q(df = 649) = 366834657.1747, p-val < .0001

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0090  0.0273  -0.3305  649  0.7411  -0.0625  0.0445    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(i2_ml(mod1),2)   
          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.44              83.05               3.64              10.76 
Code
orchard_plot(mod1, xlab = "Effect Size: Zr", group = "paperID")

4 Uni-moderator eta-regression models

Code
mod_sex <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-262.5666   525.1332   537.1332   563.9673   537.2645   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1160  0.3406    650     no         effectID 
sigma^2.2  0.0056  0.0747    197     no          paperID 
sigma^2.3  0.0173  0.1315    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 647) = 366204760.1668, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 647) = 0.7981, p-val = 0.4952

Model Results:

      estimate      se     tval   df    pval    ci.lb   ci.ub    
sexB   -0.0589  0.0425  -1.3862  647  0.1662  -0.1423  0.0245    
sexF    0.0067  0.0329   0.2037  647  0.8387  -0.0579  0.0713    
sexM    0.0013  0.0341   0.0373  647  0.9702  -0.0656  0.0682    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex)
   R2_marginal R2_conditional 
   0.004163698    0.168134011 
Code
orchard_plot(mod_sex, mod = "sex", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Code
tab_sex <- all_models(mod_sex,  mod = "sex", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_sex, here("Rdata", "tab_sex.RDS"))
Code
tab_sex <-readRDS(here("Rdata", "tab_sex.RDS"))

tab_sex
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B -0.059 -0.142 0.025 0.166 -0.795 0.678
F 0.007 -0.058 0.071 0.839 -0.728 0.741
M 0.001 -0.066 0.068 0.970 -0.734 0.736
B-F 0.066 -0.022 0.154 0.143 -0.671 0.803
B-M 0.060 -0.030 0.150 0.191 -0.677 0.798
F-M -0.005 -0.066 0.055 0.860 -0.740 0.729
Code
mod_type <- rma.mv(yi = yi,  
               V = VCV,
               mod = ~ study_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_type)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-263.3033   526.6067   536.6067   558.9761   536.7002   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1170  0.3421    650     no         effectID 
sigma^2.2  0.0045  0.0670    197     no          paperID 
sigma^2.3  0.0170  0.1303    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 648) = 366834309.7496, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 648) = 0.1628, p-val = 0.8498

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
study_typenatural        -0.0081  0.0282  -0.2864  648  0.7747  -0.0634  0.0472 
study_typesemi-natural   -0.0336  0.0609  -0.5511  648  0.5818  -0.1531  0.0860 
                          
study_typenatural         
study_typesemi-natural    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_type)
   R2_marginal R2_conditional 
  0.0004512989   0.1555024573 
Code
orchard_plot(mod_type, mod = "study_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Code
tab_type <- all_models(mod_type,  mod = "study_type", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_type, here("Rdata", "tab_type.RDS"))
Code
tab_type <-readRDS(here("Rdata", "tab_type.RDS"))

tab_type
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Natural -0.008 -0.063 0.047 0.775 -0.741 0.725
Semi-natural -0.034 -0.153 0.086 0.582 -0.774 0.707
Natural-Semi-natural -0.026 -0.143 0.092 0.671 -0.766 0.715

4.1 Higher age grouping

Code
mod_age1 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age1)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-262.1890   524.3780   536.3780   563.2120   536.5092   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1172  0.3423    650     no         effectID 
sigma^2.2  0.0045  0.0672    197     no          paperID 
sigma^2.3  0.0159  0.1263    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 647) = 364984808.3066, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 647) = 0.9677, p-val = 0.4075

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
age_class_cleanadult      -0.0043  0.0286  -0.1518  647  0.8794  -0.0606 
age_class_cleanjuvenile    0.0123  0.0472   0.2614  647  0.7939  -0.0804 
age_class_cleanmix        -0.1007  0.0626  -1.6101  647  0.1079  -0.2236 
                          ci.ub    
age_class_cleanadult     0.0519    
age_class_cleanjuvenile  0.1051    
age_class_cleanmix       0.0221    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_age1)
   R2_marginal R2_conditional 
     0.0047946      0.1528183 
Code
orchard_plot(mod_age1, mod = "age_class_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Code
tab_age1 <- all_models(mod_age1,  mod = "age_class_clean", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_age1, here("Rdata", "tab_age1.RDS"))
Code
tab_age1 <-readRDS(here("Rdata", "tab_age1.RDS"))

tab_age1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Adult -0.004 -0.061 0.052 0.879 -0.735 0.726
Juvenile 0.012 -0.080 0.105 0.794 -0.722 0.747
Mix -0.101 -0.224 0.022 0.108 -0.839 0.638
Adult-Juvenile 0.017 -0.072 0.106 0.713 -0.717 0.751
Adult-Mix -0.096 -0.217 0.025 0.118 -0.835 0.642
Juvenile-Mix -0.113 -0.254 0.027 0.115 -0.855 0.629

4.2 Finer age grouping

Code
# age_class

mod_age2 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age2)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-260.6170   521.2340   541.2340   585.8954   541.5821   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1169  0.3419    650     no         effectID 
sigma^2.2  0.0045  0.0671    197     no          paperID 
sigma^2.3  0.0170  0.1302    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 643) = 363287067.8159, p-val < .0001

Test of Moderators (coefficients 1:7):
F(df1 = 7, df2 = 643) = 0.8174, p-val = 0.5730

Model Results:

              estimate      se     tval   df    pval    ci.lb    ci.ub    
age_classA      0.0028  0.0319   0.0879  643  0.9300  -0.0597   0.0653    
age_classJ      0.0076  0.0478   0.1595  643  0.8733  -0.0863   0.1016    
age_classJA    -0.0791  0.1738  -0.4551  643  0.6492  -0.4203   0.2622    
age_classJY    -0.2024  0.0945  -2.1412  643  0.0326  -0.3881  -0.0168  * 
age_classJYA   -0.0223  0.0889  -0.2504  643  0.8023  -0.1968   0.1523    
age_classY     -0.0267  0.0616  -0.4337  643  0.6647  -0.1478   0.0943    
age_classYA    -0.0305  0.0460  -0.6615  643  0.5085  -0.1209   0.0600    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_age2)
   R2_marginal R2_conditional 
    0.01008384     0.16366782 
Code
orchard_plot(mod_age2, mod = "age_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Code
# table not created for this (not meaningful for some levels)
Code
mod_gen <- rma.mv(yi = yi, 
                V = VCV,
               mod = ~ whose_fitness - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_gen)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-262.7747   525.5494   535.5494   557.9188   535.6428   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1163  0.3411    650     no         effectID 
sigma^2.2  0.0030  0.0545    197     no          paperID 
sigma^2.3  0.0200  0.1414    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 648) = 366295125.0150, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 648) = 0.9896, p-val = 0.3723

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
whose_fitnessindividual   -0.0208  0.0286  -0.7283  648  0.4667  -0.0770 
whose_fitnessoffspring     0.0333  0.0446   0.7479  648  0.4548  -0.0542 
                          ci.ub    
whose_fitnessindividual  0.0353    
whose_fitnessoffspring   0.1208    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
orchard_plot(mod_gen, mod = "whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Code
tab_gen <- all_models(mod_gen,  mod = "whose_fitness", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_gen, here("Rdata", "tab_gen.RDS"))
Code
tab_gen <-readRDS(here("Rdata", "tab_gen.RDS"))

tab_gen
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Individual -0.021 -0.077 0.035 0.467 -0.756 0.714
Offspring 0.033 -0.054 0.121 0.455 -0.705 0.771
Individual-Offspring 0.054 -0.026 0.134 0.184 -0.683 0.791

4.3 Higher level

Code
mod_fit1 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_higher_level - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit1)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-262.7944   525.5889   537.5889   564.4230   537.7201   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1172  0.3424    650     no         effectID 
sigma^2.2  0.0053  0.0730    197     no          paperID 
sigma^2.3  0.0152  0.1232    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 647) = 365633577.2030, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 647) = 0.6085, p-val = 0.6097

Model Results:

                                      estimate      se     tval   df    pval 
fitness_higher_levelindirect fitness   -0.2262  0.1684  -1.3435  647  0.1796 
fitness_higher_levelreproduction       -0.0047  0.0308  -0.1529  647  0.8785 
fitness_higher_levelsurvival           -0.0075  0.0319  -0.2359  647  0.8136 
                                        ci.lb   ci.ub    
fitness_higher_levelindirect fitness  -0.5568  0.1044    
fitness_higher_levelreproduction      -0.0651  0.0557    
fitness_higher_levelsurvival          -0.0703  0.0552    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_fit1)
   R2_marginal R2_conditional 
   0.002697139    0.151162291 
Code
orchard_plot(mod_fit1, mod = "fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
tab_fit1 <- all_models(mod_fit1,  mod = "fitness_higher_level", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_fit1, here("Rdata", "tab_fit1.RDS"))
Code
tab_fit1 <-readRDS(here("Rdata", "tab_fit1.RDS"))

tab_fit1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Indirect fitness -0.226 -0.557 0.104 0.180 -1.026 0.574
Reproduction -0.005 -0.065 0.056 0.878 -0.736 0.727
Survival -0.008 -0.070 0.055 0.814 -0.739 0.724
Indirect fitness-Reproduction 0.222 -0.111 0.554 0.191 -0.579 1.022
Indirect fitness-Survival 0.219 -0.113 0.551 0.197 -0.582 1.020
Reproduction-Survival -0.003 -0.063 0.057 0.926 -0.734 0.728

4.4 Finer level

  • need to think what to do with classes less than 5 studies
Code
mod_fit2 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_metric_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit2)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-261.8909   523.7818   545.7818   594.8923   546.2009   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1170  0.3420    650     no         effectID 
sigma^2.2  0.0035  0.0594    197     no          paperID 
sigma^2.3  0.0193  0.1388    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 642) = 364423843.9729, p-val < .0001

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 642) = 0.5448, p-val = 0.8229

Model Results:

                                                           estimate      se 
fitness_metric_cleanage at first reproduction               -0.0066  0.1123 
fitness_metric_cleanindirect fitness                        -0.2244  0.1694 
fitness_metric_cleanlifetime breeding success                0.0032  0.0368 
fitness_metric_cleanlifetime reproductive success           -0.0209  0.0410 
fitness_metric_cleanoffspring reproduction                   0.0445  0.1098 
fitness_metric_cleanoffspring survival                       0.0321  0.0466 
fitness_metric_cleanreproductive lifespan and/or attempts   -0.0752  0.1402 
fitness_metric_cleansurvival                                -0.0301  0.0361 
                                                              tval   df    pval 
fitness_metric_cleanage at first reproduction              -0.0588  642  0.9531 
fitness_metric_cleanindirect fitness                       -1.3248  642  0.1857 
fitness_metric_cleanlifetime breeding success               0.0883  642  0.9296 
fitness_metric_cleanlifetime reproductive success          -0.5104  642  0.6100 
fitness_metric_cleanoffspring reproduction                  0.4054  642  0.6853 
fitness_metric_cleanoffspring survival                      0.6893  642  0.4909 
fitness_metric_cleanreproductive lifespan and/or attempts  -0.5366  642  0.5918 
fitness_metric_cleansurvival                               -0.8348  642  0.4041 
                                                             ci.lb   ci.ub    
fitness_metric_cleanage at first reproduction              -0.2271  0.2139    
fitness_metric_cleanindirect fitness                       -0.5571  0.1082    
fitness_metric_cleanlifetime breeding success              -0.0690  0.0755    
fitness_metric_cleanlifetime reproductive success          -0.1014  0.0596    
fitness_metric_cleanoffspring reproduction                 -0.1710  0.2600    
fitness_metric_cleanoffspring survival                     -0.0594  0.1237    
fitness_metric_cleanreproductive lifespan and/or attempts  -0.3505  0.2001    
fitness_metric_cleansurvival                               -0.1010  0.0407    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_fit2)
   R2_marginal R2_conditional 
   0.006544495    0.168509539 
Code
orchard_plot(mod_fit2, mod = "fitness_metric_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
# table not created for this (not meaningful for some levels)
Code
mod_disp1 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ dispersal_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_disp1)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-263.0695   526.1391   538.1391   564.9732   538.2703   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1170  0.3421    650     no         effectID 
sigma^2.2  0.0041  0.0643    197     no          paperID 
sigma^2.3  0.0176  0.1327    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 647) = 366209372.4133, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 647) = 0.3199, p-val = 0.8110

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
dispersal_typeB          -0.0420  0.0805  -0.5217  647  0.6020  -0.2000  0.1160 
dispersal_typebreeding    0.0141  0.0405   0.3493  647  0.7270  -0.0653  0.0936 
dispersal_typenatal      -0.0173  0.0299  -0.5784  647  0.5632  -0.0760  0.0414 
                          
dispersal_typeB           
dispersal_typebreeding    
dispersal_typenatal       

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_disp1)
   R2_marginal R2_conditional 
   0.001665491    0.158129749 
Code
orchard_plot(mod_disp1, mod = "dispersal_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
tab_disp1 <- all_models(mod_disp1,  mod = "dispersal_type", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_disp1, here("Rdata", "tab_disp1.RDS"))
Code
tab_disp1 <-readRDS(here("Rdata", "tab_disp1.RDS"))

tab_disp1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B -0.042 -0.200 0.116 0.602 -0.790 0.706
Breeding 0.014 -0.065 0.094 0.727 -0.722 0.750
Natal -0.017 -0.076 0.041 0.563 -0.751 0.717
B-Breeding 0.056 -0.111 0.223 0.510 -0.694 0.806
B-Natal 0.025 -0.134 0.183 0.760 -0.724 0.773
Breeding-Natal -0.031 -0.109 0.046 0.427 -0.767 0.704
Code
mod_disp2 <- rma.mv(yi = yi, V = vi,
               mod = ~ dispersal_phase - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_disp2)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-208.9178   417.8357   427.8357   450.2051   427.9291   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0681  0.2610    650     no         effectID 
sigma^2.2  0.0051  0.0713    197     no          paperID 
sigma^2.3  0.0169  0.1299    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 648) = 7264.6862, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 648) = 1.1574, p-val = 0.3150

Model Results:

                           estimate      se     tval   df    pval    ci.lb 
dispersal_phasedispersal    -0.0475  0.0340  -1.3987  648  0.1624  -0.1142 
dispersal_phasesettlement    0.0084  0.0216   0.3876  648  0.6984  -0.0341 
                            ci.ub    
dispersal_phasedispersal   0.0192    
dispersal_phasesettlement  0.0508    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_disp2)
   R2_marginal R2_conditional 
   0.005618405    0.248040842 
Code
orchard_plot(mod_disp2, mod = "dispersal_phase", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
tab_disp2 <- all_models(mod_disp2,  mod = "dispersal_phase", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_disp2, here("Rdata", "tab_disp2.RDS"))
Code
tab_disp2 <-readRDS(here("Rdata", "tab_disp2.RDS"))

tab_disp2
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Dispersal -0.047 -0.114 0.019 0.162 -0.641 0.546
Settlement 0.008 -0.034 0.051 0.698 -0.582 0.599
Dispersal-Settlement 0.048 -0.033 0.129 0.244 -0.689 0.785
Code
mod_class <- rma.mv(yi = yi, V = vi,
               mod = ~ species_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_class)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-207.9232   415.8463   433.8463   474.0556   434.1302   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0687  0.2621    650     no         effectID 
sigma^2.2  0.0054  0.0735    197     no          paperID 
sigma^2.3  0.0158  0.1257    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 644) = 7073.8786, p-val < .0001

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 644) = 0.5936, p-val = 0.7356

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
species_classActinopterygii    0.0346  0.1398   0.2475  644  0.8046  -0.2400 
species_classArachnida         0.0464  0.2227   0.2085  644  0.8349  -0.3909 
species_classAves             -0.0160  0.0260  -0.6133  644  0.5399  -0.0671 
species_classInsecta           0.1353  0.1216   1.1119  644  0.2666  -0.1036 
species_classMammalia         -0.0143  0.0312  -0.4573  644  0.6476  -0.0756 
species_classReptilia          0.1466  0.1147   1.2788  644  0.2014  -0.0785 
                              ci.ub    
species_classActinopterygii  0.3092    
species_classArachnida       0.4838    
species_classAves            0.0352    
species_classInsecta         0.3741    
species_classMammalia        0.0470    
species_classReptilia        0.3718    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_class)
   R2_marginal R2_conditional 
    0.01183027     0.24481984 
Code
orchard_plot(mod_class, mod = "species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
tab_class <- all_models(mod_class,  mod = "species_class", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_class, here("Rdata", "tab_class.RDS"))
Code
tab_class <-readRDS(here("Rdata", "tab_class.RDS"))

tab_class
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Actinopterygii 0.035 -0.240 0.309 0.805 -0.615 0.684
Arachnida 0.046 -0.391 0.484 0.835 -0.687 0.780
Aves -0.016 -0.067 0.035 0.540 -0.607 0.575
Insecta 0.135 -0.104 0.374 0.267 -0.500 0.771
Mammalia -0.014 -0.076 0.047 0.648 -0.606 0.578
Reptilia 0.147 -0.079 0.372 0.201 -0.484 0.777
Actinopterygii-Arachnida 0.048 -0.535 0.631 0.871 -0.893 0.990
Actinopterygii-Aves -0.017 -0.326 0.291 0.912 -0.818 0.784
Actinopterygii-Insecta 0.116 -0.292 0.523 0.577 -0.728 0.960
Actinopterygii-Mammalia -0.032 -0.342 0.279 0.842 -0.833 0.770
Actinopterygii-Reptilia 0.084 -0.299 0.468 0.667 -0.749 0.917
Arachnida-Aves -0.065 -0.568 0.437 0.798 -0.959 0.828
Arachnida-Insecta 0.068 -0.501 0.636 0.815 -0.865 1.000
Arachnida-Mammalia -0.080 -0.583 0.424 0.756 -0.974 0.815
Arachnida-Reptilia 0.036 -0.516 0.588 0.898 -0.886 0.958
Aves-Insecta 0.133 -0.146 0.413 0.349 -0.657 0.923
Aves-Mammalia -0.014 -0.105 0.077 0.759 -0.759 0.730
Aves-Reptilia 0.101 -0.143 0.345 0.415 -0.677 0.880
Insecta-Mammalia -0.147 -0.429 0.134 0.304 -0.938 0.643
Insecta-Reptilia -0.032 -0.393 0.329 0.862 -0.854 0.791
Mammalia-Reptilia 0.116 -0.129 0.360 0.355 -0.663 0.894
Code
mod_design <- rma.mv(yi = yi, V = vi,
               mod = ~ study_design - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_design)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-209.1365   418.2730   428.2730   450.6425   428.3665   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0692  0.2631    650     no         effectID 
sigma^2.2  0.0080  0.0892    197     no          paperID 
sigma^2.3  0.0103  0.1016    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 648) = 7226.9352, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 648) = 0.7488, p-val = 0.4733

Model Results:

                      estimate      se     tval   df    pval    ci.lb   ci.ub 
study_designcontrast    0.0069  0.0203   0.3390  648  0.7348  -0.0330  0.0468 
study_designgradient   -0.0390  0.0343  -1.1353  648  0.2567  -0.1064  0.0284 
                        
study_designcontrast    
study_designgradient    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_design)
   R2_marginal R2_conditional 
   0.004390445    0.212426304 
Code
orchard_plot(mod_design, mod = "study_design", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
tab_design <- all_models(mod_design,  mod = "study_design", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_design, here("Rdata", "tab_design.RDS"))
Code
tab_design <-readRDS(here("Rdata", "tab_design.RDS"))

tab_design
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Contrast 0.007 -0.033 0.047 0.735 -0.575 0.589
Gradient -0.039 -0.106 0.028 0.257 -0.624 0.546
Contrast-Gradient -0.025 -0.107 0.057 0.553 -0.758 0.709

4.5 Normal analysis

Code
mod_focus <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_focus)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-263.3982   526.7963   536.7963   559.1658   536.8898   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1169  0.3419    650     no         effectID 
sigma^2.2  0.0053  0.0731    197     no          paperID 
sigma^2.3  0.0161  0.1269    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 648) = 366616051.8753, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 648) = 0.1217, p-val = 0.8855

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub    
fitness_main_focusN    0.0007  0.0424   0.0159  648  0.9873  -0.0825  0.0838    
fitness_main_focusY   -0.0129  0.0288  -0.4472  648  0.6549  -0.0695  0.0437    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_focus)
   R2_marginal R2_conditional 
  0.0002449745   0.1551757994 
Code
orchard_plot(mod_focus, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Code
tab_focus <- all_models(mod_focus,  mod = "fitness_main_focus", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_focus, here("Rdata", "tab_focus.RDS"))
Code
tab_focus <-readRDS(here("Rdata", "tab_focus.RDS"))

tab_focus
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
N 0.001 -0.082 0.084 0.987 -0.734 0.736
Y -0.013 -0.069 0.044 0.655 -0.745 0.720
N-Y -0.014 -0.094 0.067 0.740 -0.748 0.721

4.6 Absolute value analysis

Code
# homo
mod_focusA <- rma.mv(yi = abs_yi2, 
                V = VCV_abs,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_focusA)

Multivariate Meta-Analysis Model (k = 650; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
-18.6194   37.2388   47.2388   69.6083   47.3323   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0529  0.2300    650     no         effectID 
sigma^2.2  0.0112  0.1060    197     no          paperID 
sigma^2.3  0.0000  0.0000    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 648) = 328986653.1626, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 648) = 2.7837, p-val = 0.0626

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub    
fitness_main_focusN   -0.0042  0.0299  -0.1392  648  0.8894  -0.0630  0.0546    
fitness_main_focusY    0.0419  0.0194   2.1641  648  0.0308   0.0039  0.0800  * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_focusA)
   R2_marginal R2_conditional 
   0.006071876    0.180198901 
Code
# hetero
mod_focusA1 <- rma.mv(yi = abs_yi2, 
                 V = VCV_abs,
                 mod = ~ fitness_main_focus,
                 data = dat, 
                 random = list(
                   ~ fitness_main_focus | effectID,
                   ~ 1 | paperID,
                   ~ 1 | species_cleaned),
                 struct = "DIAG",
                 method = "ML",
                 test = "t",
                 sparse = TRUE,
                 control = list(optimizer = "Nelder-Mead"))
summary(mod_focusA1)

Multivariate Meta-Analysis Model (k = 650; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
  -7.7128  9099.6507    27.4257    54.2875    27.5563   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0099  0.0995    197     no          paperID 
sigma^2.2  0.0000  0.0002    140     no  species_cleaned 

outer factor: effectID           (nlvls = 650)
inner factor: fitness_main_focus (nlvls = 2)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.0294  0.1716    158     no      N 
tau^2.2    0.0613  0.2475    492     no      Y 

Test for Residual Heterogeneity:
QE(df = 648) = 328986653.1626, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 648) = 2.4929, p-val = 0.1148

Model Results:

                     estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt                0.0146  0.0256  0.5694  648  0.5693  -0.0357  0.0648    
fitness_main_focusY    0.0412  0.0261  1.5789  648  0.1148  -0.0100  0.0924    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
#r2_ml(mod_focusA1)

# homo
orchard_plot(mod_focusA, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Code
# hetero
orchard_plot(mod_focusA1, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Code
tab_focusA1 <- all_models(mod_focusA1,  mod = "fitness_main_focus", type = "hetero")

# saving tab_sex as RDS
saveRDS(tab_focus, here("Rdata", "tab_focusA1.RDS"))
Code
tab_focusA1 <-readRDS(here("Rdata", "tab_focusA1.RDS"))

tab_focusA1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
N 0.001 -0.082 0.084 0.987 -0.734 0.736
Y -0.013 -0.069 0.044 0.655 -0.745 0.720
N-Y -0.014 -0.094 0.067 0.740 -0.748 0.721

5 Multi-moderator-regression model

Code
# use "ML" for model selection
#######################
# Mulit-variable models
#######################

# we need - to
mod_full <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1 + sex +
                age_class_clean + whose_fitness + 
                fitness_higher_level +
                dispersal_type + dispersal_phase,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              method = "ML",
              test = "t",
              sparse = TRUE)
summary(mod_full)

Multivariate Meta-Analysis Model (k = 650; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-259.7157  9198.5103   547.4314   610.1090   548.0928   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1145  0.3383    650     no         effectID 
sigma^2.2  0.0007  0.0270    197     no          paperID 
sigma^2.3  0.0234  0.1531    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 639) = 361692892.8223, p-val < .0001

Test of Moderators (coefficients 2:11):
F(df1 = 10, df2 = 639) = 0.9415, p-val = 0.4941

Model Results:

                                  estimate      se     tval   df    pval 
intrcpt                            -0.2864  0.1860  -1.5402  639  0.1240 
sexF                                0.0500  0.0485   1.0312  639  0.3028 
sexM                                0.0461  0.0494   0.9318  639  0.3518 
age_class_cleanjuvenile             0.0577  0.0529   1.0915  639  0.2755 
age_class_cleanmix                 -0.0635  0.0661  -0.9604  639  0.3372 
whose_fitnessoffspring              0.0494  0.0489   1.0102  639  0.3128 
fitness_higher_levelreproduction    0.2197  0.1705   1.2885  639  0.1980 
fitness_higher_levelsurvival        0.2167  0.1731   1.2518  639  0.2111 
dispersal_typebreeding              0.0049  0.0873   0.0566  639  0.9548 
dispersal_typenatal                -0.0294  0.0821  -0.3582  639  0.7203 
dispersal_phasesettlement           0.0363  0.0511   0.7113  639  0.4772 
                                    ci.lb   ci.ub    
intrcpt                           -0.6516  0.0788    
sexF                              -0.0452  0.1453    
sexM                              -0.0510  0.1432    
age_class_cleanjuvenile           -0.0461  0.1615    
age_class_cleanmix                -0.1934  0.0664    
whose_fitnessoffspring            -0.0466  0.1454    
fitness_higher_levelreproduction  -0.1151  0.5546    
fitness_higher_levelsurvival      -0.1232  0.5567    
dispersal_typebreeding            -0.1665  0.1764    
dispersal_typenatal               -0.1907  0.1319    
dispersal_phasesettlement         -0.0640  0.1367    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_full)*100, 2)
   R2_marginal R2_conditional 
          1.62          18.77 
Code
# multi-model selection
candidates <- dredge(mod_full, trace = 2)

# saving tab_sex as RDS
saveRDS(candidates, here("Rdata", "candidates.RDS"))
Code
candidates <-readRDS(here("Rdata", "candidates.RDS"))

candidates
Global model call: rma.mv(yi = yi, V = VCV, mods = ~1 + sex + age_class_clean + 
    whose_fitness + fitness_higher_level + dispersal_type + dispersal_phase, 
    random = list(~1 | effectID, ~1 | paperID, ~1 | species_cleaned), 
    data = dat, method = "ML", test = "t", sparse = TRUE)
---
Model selection table 
   (Int) age_cls_cln dsp_phs dsp_typ ftn_hgh_lvl sex whs_ftn df   logLik  AICc
1      +                                                      4 -264.109 536.3
33     +                                                   +  5 -263.311 536.7
3      +                   +                                  5 -263.495 537.1
2      +           +                                          6 -262.717 537.6
35     +                   +                               +  6 -262.796 537.7
34     +           +                                       +  7 -262.007 538.2
17     +                                           +          6 -263.045 538.2
4      +           +       +                                  7 -262.071 538.3
9      +                                       +              6 -263.248 538.6
41     +                                       +           +  7 -262.356 538.9
49     +                                           +       +  7 -262.438 539.1
36     +           +       +                               +  8 -261.430 539.1
11     +                   +                   +              7 -262.513 539.2
5      +                           +                          6 -263.726 539.6
37     +                           +                       +  7 -262.794 539.8
18     +           +                               +          8 -261.792 539.8
19     +                   +                       +          7 -262.821 539.8
10     +           +                           +              8 -261.821 539.9
43     +                   +                   +           +  8 -261.988 540.2
25     +                                       +   +          8 -262.084 540.4
42     +           +                           +           +  9 -261.057 540.4
12     +           +       +                   +              9 -261.067 540.4
51     +                   +                       +       +  8 -262.240 540.7
7      +                   +       +                          7 -263.282 540.7
50     +           +                               +       +  9 -261.232 540.7
6      +           +               +                          8 -262.446 541.1
39     +                   +       +                       +  8 -262.456 541.1
20     +           +       +                       +          9 -261.439 541.2
38     +           +               +                       +  9 -261.543 541.4
57     +                                       +   +       +  9 -261.582 541.4
44     +           +       +                   +           + 10 -260.597 541.5
27     +                   +                   +   +          9 -261.740 541.8
21     +                           +               +          8 -262.792 541.8
26     +           +                           +   +         10 -260.822 542.0
8      +           +       +       +                          9 -261.876 542.0
13     +                           +           +              8 -262.947 542.1
52     +           +       +                       +       + 10 -260.907 542.2
45     +                           +           +           +  9 -261.981 542.2
53     +                           +               +       +  9 -262.077 542.4
40     +           +       +       +                       + 10 -261.051 542.4
15     +                   +       +           +              9 -262.341 543.0
58     +           +                           +   +       + 11 -260.343 543.1
28     +           +       +                   +   +         11 -260.363 543.1
59     +                   +                   +   +       + 10 -261.401 543.1
22     +           +               +               +         10 -261.537 543.4
14     +           +               +           +             10 -261.591 543.5
23     +                   +       +               +          9 -262.646 543.6
46     +           +               +           +           + 11 -260.649 543.7
47     +                   +       +           +           + 10 -261.719 543.8
54     +           +               +               +       + 11 -260.816 544.0
29     +                           +           +   +         10 -261.866 544.1
16     +           +       +       +           +             11 -260.880 544.2
55     +                   +       +               +       + 10 -261.961 544.3
60     +           +       +                   +   +       + 12 -260.047 544.6
24     +           +       +       +               +         11 -261.238 544.9
48     +           +       +       +           +           + 12 -260.249 545.0
61     +                           +           +   +       + 11 -261.302 545.0
56     +           +       +       +               +       + 12 -260.543 545.6
31     +                   +       +           +   +         11 -261.594 545.6
30     +           +               +           +   +         12 -260.584 545.7
62     +           +               +           +   +       + 13 -259.966 546.5
63     +                   +       +           +   +       + 12 -261.181 546.9
32     +           +       +       +           +   +         13 -260.161 546.9
64     +           +       +       +           +   +       + 14 -259.716 548.1
   delta weight
1   0.00  0.110
33  0.43  0.089
3   0.80  0.074
2   1.28  0.058
35  1.44  0.054
34  1.91  0.042
17  1.94  0.042
4   2.03  0.040
9   2.35  0.034
41  2.61  0.030
49  2.77  0.028
36  2.80  0.027
11  2.92  0.026
5   3.30  0.021
37  3.48  0.019
18  3.53  0.019
19  3.54  0.019
10  3.59  0.018
43  3.92  0.016
25  4.11  0.014
42  4.12  0.014
12  4.13  0.014
51  4.42  0.012
7   4.46  0.012
50  4.47  0.012
6   4.84  0.010
39  4.86  0.010
20  4.88  0.010
38  5.09  0.009
57  5.17  0.008
44  5.26  0.008
27  5.48  0.007
21  5.53  0.007
26  5.71  0.006
8   5.75  0.006
13  5.84  0.006
52  5.88  0.006
45  5.96  0.006
53  6.15  0.005
40  6.16  0.005
15  6.68  0.004
58  6.82  0.004
28  6.86  0.004
59  6.86  0.004
22  7.14  0.003
14  7.25  0.003
23  7.29  0.003
46  7.43  0.003
47  7.50  0.003
54  7.77  0.002
29  7.79  0.002
16  7.89  0.002
55  7.98  0.002
60  8.30  0.002
24  8.61  0.001
48  8.71  0.001
61  8.74  0.001
56  9.30  0.001
31  9.32  0.001
30  9.38  0.001
62 10.22  0.001
63 10.57  0.001
32 10.61  0.001
64 11.81  0.000
Models ranked by AICc(x) 
Code
# displays delta AICc <2
candidates_aic2 <- subset(candidates, delta < 2) 
# model averaging
mr_averaged_aic2 <- summary(model.avg(candidates, delta < 2)) 

mr_averaged_aic2

Call:
model.avg(object = candidates, subset = delta < 2)

Component model call: 
rma.mv(yi = yi, V = VCV, mods = ~<7 unique rhs>, random = list(~1 | 
     effectID, ~1 | paperID, ~1 | species_cleaned), data = dat, method = ML, 
     test = t, sparse = TRUE)

Component models: 
       df  logLik   AICc delta weight
(Null)  4 -264.11 536.28  0.00   0.24
4       5 -263.31 536.71  0.43   0.19
2       5 -263.49 537.08  0.80   0.16
1       6 -262.72 537.57  1.28   0.12
24      6 -262.80 537.72  1.44   0.11
14      7 -262.01 538.19  1.91   0.09
3       6 -263.05 538.22  1.94   0.09

Term codes: 
age_class_clean dispersal_phase             sex   whose_fitness 
              1               2               3               4 

Model-averaged coefficients:  
(full average) 
                           Estimate Std. Error z value Pr(>|z|)
intrcpt                   -0.026596   0.039875   0.667    0.505
whose_fitnessoffspring     0.020470   0.036076   0.567    0.570
dispersal_phasesettlement  0.012379   0.029609   0.418    0.676
age_class_cleanjuvenile    0.004109   0.022391   0.183    0.854
age_class_cleanmix        -0.020044   0.047793   0.419    0.675
sexF                       0.005767   0.022722   0.254    0.800
sexM                       0.005266   0.021666   0.243    0.808
 
(conditional average) 
                          Estimate Std. Error z value Pr(>|z|)
intrcpt                   -0.02660    0.03988   0.667    0.505
whose_fitnessoffspring     0.05191    0.04085   1.271    0.204
dispersal_phasesettlement  0.04553    0.04141   1.099    0.272
age_class_cleanjuvenile    0.01917    0.04528   0.423    0.672
age_class_cleanmix        -0.09351    0.06154   1.520    0.129
sexF                       0.06468    0.04450   1.453    0.146
sexM                       0.05906    0.04569   1.292    0.196
Code
# relative importance of each predictor for all the models
importance <- sw(candidates)

importance
                     whose_fitness dispersal_phase age_class_clean
Sum of weights:      0.42          0.37            0.33           
N containing models:   32            32              32           
                     fitness_higher_level sex  dispersal_type
Sum of weights:      0.24                 0.23 0.15          
N containing models:   32                   32   32          

6 Publication bias

Code
# write conditional statement here

# if each group n is avaiable - assume n/2
dat$effectN <- ifelse(is.na(dat$n_group_1), (dat$n/2)*2/dat$n,  
                      (dat$n_group_1 * dat$n_group_2) / (dat$n_group_1 + dat$n_group_2))
  
dat$sqeffectN <- sqrt(dat$effectN)

mod_egger <- rma.mv(yi = yi, 
                    V = VCV,
                mod = ~ sqeffectN,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                # ~ 1 | phylogeny),
                # R= list(phylogeny = cor_tree),
                test = "t",
                sparse = TRUE)
summary(mod_egger)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-263.4161   526.8323   536.8323   559.2017   536.9257   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1174  0.3426    650     no         effectID 
sigma^2.2  0.0055  0.0743    197     no          paperID 
sigma^2.3  0.0149  0.1223    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 648) = 366820515.9753, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 648) = 0.0004, p-val = 0.9839

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -0.0095  0.0381  -0.2485  648  0.8038  -0.0842  0.0653    
sqeffectN    0.0001  0.0044   0.0202  648  0.9839  -0.0085  0.0087    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_egger)*100, 2)
   R2_marginal R2_conditional 
          0.00          14.86 
Code
small <- bubble_plot(mod_egger,
                     mod = "sqeffectN",
                     group = "paperID",
                     xlab = "sqrt(Effective N)",
                     ylab = "Effect Size: Zr",
                     g = TRUE)

small

Code
# creating publication year from "reference"
dat$reference <- as.character(dat$reference)
dat$year <- with(dat, substr(reference, nchar(reference)-4, nchar(reference)))
dat$year <- as.integer(ifelse(dat$year == "2017a" | dat$year == "2017b", 2017, dat$year))
# decline effect
mod_dec <- rma.mv(yi = yi, V = vi,
                     mod = ~ year,
                     data = dat, 
                     random = list(
                       ~ 1 | effectID,
                       ~ 1 | paperID,
                       ~ 1 | species_cleaned),
                     # ~ 1 | phylogeny),
                     # R= list(phylogeny = cor_tree),
                     test = "t",
                     sparse = TRUE)
summary(mod_dec)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-208.8588   417.7177   427.7177   450.0871   427.8111   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0689  0.2624    650     no         effectID 
sigma^2.2  0.0056  0.0750    197     no          paperID 
sigma^2.3  0.0141  0.1189    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 648) = 7252.6587, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 648) = 1.9742, p-val = 0.1605

Model Results:

         estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt   -3.9042  2.7751  -1.4069  648  0.1599  -9.3535  1.5451    
year       0.0019  0.0014   1.4051  648  0.1605  -0.0008  0.0047    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_dec)*100, 2)
   R2_marginal R2_conditional 
          0.62          22.77 
Code
decline <- bubble_plot(mod_dec,
                       mod = "year",
                       group = "paperID",
                       xlab = "Publication year",
                       ylab = "Effect Size: Zr",
                       g = TRUE)
decline 

Code
# All together

mod_comb <- rma.mv(yi = yi, 
                   V = VCV,
                   mod = ~ year + sqeffectN,
                   data = dat, 
                   random = list(
                     ~ 1 | effectID,
                     ~ 1 | paperID,
                     ~ 1 | species_cleaned),
                   # ~ 1 | phylogeny),
                   # R= list(phylogeny = cor_tree),
                   test = "t",
                   sparse = TRUE)
                   
summary(mod_comb)

Multivariate Meta-Analysis Model (k = 650; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-262.7944   525.5889   537.5889   564.4229   537.7201   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1166  0.3415    650     no         effectID 
sigma^2.2  0.0053  0.0726    197     no          paperID 
sigma^2.3  0.0169  0.1302    140     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 647) = 366814490.8912, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 647) = 0.4654, p-val = 0.6281

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -3.0343  3.1339  -0.9682  647  0.3333  -9.1881  3.1196    
year         0.0015  0.0016   0.9647  647  0.3350  -0.0016  0.0046    
sqeffectN   -0.0006  0.0044  -0.1292  647  0.8973  -0.0093  0.0081    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_comb)*100, 2)
   R2_marginal R2_conditional 
          0.23          16.20 
Code
# small-study
bubble_plot(mod_egger,
            mod = "sqeffectN",
            group = "paperID",
            xlab = "sqrt(Effective N)",
            ylab = "Effect Size: Zr",
            g = TRUE)

Code
# decline
bubble_plot(mod_comb,
            mod = "year",
            group = "paperID",
            xlab = "Publication year",
            ylab = "Effect Size: Zr",
            g = TRUE)

Code
dat <- dat %>%
  mutate(leave_out = reference)

dat$leave_out <- as.factor(dat$leave_out)


LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$leave_out))) {
  temp_dat <- dat %>%
    filter(leave_out != levels(dat$leave_out)[i])
  
  VCV_leaveout <- vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, rho = 0.5)
  
  LeaveOneOut_effectsize[[i]] <-  rma.mv(yi = yi,
                                         V = VCV_leaveout, 
                                         random = list(
                                           ~ 1 | effectID,
                                           ~ 1 | paperID,
                                           ~ 1 | species_cleaned),
                                         # ~ 1 | phylogeny),
                                         # R= list(phylogeny = cor_tree),
                                         test = "t",
                                         sparse = TRUE,
                                         data = temp_dat[temp_dat$leave_out != levels(temp_dat$leave_out)[i], ])
}

# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(model) {
  df <- data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub)
  return(df)
}

# using dplyr to form data frame
MA_oneout <- lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>%
  bind_rows %>%
  mutate(left_out = levels(dat$leave_out))


# telling ggplot to stop reordering factors
MA_oneout$left_out <- as.factor(MA_oneout$left_out)
MA_oneout$left_out <- factor(MA_oneout$left_out, levels = MA_oneout$left_out)

# saving the runs
saveRDS(MA_oneout, here("Rdata", "MA_oneout.RDS"))
Code
MA_oneout <- readRDS(here("Rdata", "MA_oneout.RDS"))

# plotting
leaveoneout <- ggplot(MA_oneout) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +

  geom_hline(yintercept = mod1$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
  geom_hline(yintercept = mod1$b, lty = 1, lwd = 0.75, colour = "black") + 
  geom_hline(yintercept = mod1$ci.ub,
             lty = 3, lwd = 0.75, colour = "black") + 
  geom_pointrange(aes(x = left_out, y = est,
                      ymin = lower, ymax = upper)) + 
  xlab("Study left out") + 
  ylab("Zr (effect size), 95% CI") +
  coord_flip() + 
  theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))

leaveoneout

8 R Session Informtion

Code
# pander for making it look nicer
sessionInfo() %>% pander()

R version 4.3.1 (2023-06-16)

Platform: aarch64-apple-darwin20 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: orchaRd(v.2.0), here(v.1.0.1), patchwork(v.1.1.3), ape(v.5.7-1), metafor(v.4.2-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), Matrix(v.1.6-1), pander(v.0.6.5), gridExtra(v.2.3), emmeans(v.1.8.8), MuMIn(v.1.47.5), kableExtra(v.1.3.4), lubridate(v.1.9.2), forcats(v.1.0.0), stringr(v.1.5.0), dplyr(v.1.1.3), purrr(v.1.0.2), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.4.4) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): tidyselect(v.1.2.0), viridisLite(v.0.4.2), farver(v.2.1.1), vipor(v.0.4.5), latex2exp(v.0.9.6), fastmap(v.1.1.1), TH.data(v.1.1-2), pacman(v.0.5.1), mathjaxr(v.1.6-0), digest(v.0.6.34), estimability(v.1.4.1), timechange(v.0.2.0), lifecycle(v.1.0.4), survival(v.3.5-7), magrittr(v.2.0.3), compiler(v.4.3.1), rlang(v.1.1.3), tools(v.4.3.1), utf8(v.1.2.3), yaml(v.2.3.8), knitr(v.1.43), labeling(v.0.4.3), htmlwidgets(v.1.6.2), xml2(v.1.3.5), multcomp(v.1.4-25), withr(v.2.5.2), grid(v.4.3.1), stats4(v.4.3.1), fansi(v.1.0.4), xtable(v.1.8-4), colorspace(v.2.1-0), scales(v.1.2.1), MASS(v.7.3-60), cli(v.3.6.2), mvtnorm(v.1.2-3), rmarkdown(v.2.24), generics(v.0.1.3), rstudioapi(v.0.15.0), httr(v.1.4.7), tzdb(v.0.4.0), ggbeeswarm(v.0.7.2), splines(v.4.3.1), rvest(v.1.0.3), parallel(v.4.3.1), vctrs(v.0.6.5), webshot(v.0.5.5), sandwich(v.3.0-2), jsonlite(v.1.8.8), hms(v.1.1.3), beeswarm(v.0.4.0), systemfonts(v.1.0.4), glue(v.1.7.0), codetools(v.0.2-19), stringi(v.1.7.12), gtable(v.0.3.4), munsell(v.0.5.0), pillar(v.1.9.0), htmltools(v.0.5.7), R6(v.2.5.1), rprojroot(v.2.0.3), evaluate(v.0.21), lattice(v.0.21-8), Rcpp(v.1.0.12), svglite(v.2.1.1), coda(v.0.19-4), nlme(v.3.1-163), mgcv(v.1.9-0), xfun(v.0.40), zoo(v.1.8-12) and pkgconfig(v.2.0.3)